FROM HOPF TO NEIMARK SACKER BIFURCATION: 
A COMPUTATIONAL ALGORITHM 



GERALD MOORE* 

Abstract . We construct an algorithm for approximating the invariant tori created at a Neimark— 
Sacker bifurcation point. It is based on the same philosophy as many algorithms for approximating 
the periodic orbits created at a Hopf bifurcation point, i.e. a Fourier spectral method. For Neimark— 
Sacker bifurcation, however, we use a simple parametrisation of the tori in order to determine low- 
order approximations, and then utilise the information contained therein to develop a more general 
parametrisation suitable for computing higher-order approximations. Different algorithms, applicable 
to either autonomous or periodically-forced systems of differential equations, are obtained. 
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1. Introduction. In this paper we consider both nonhnear autonomous systems 

— =F{x,X) F:M"xMk^M", (1.1) 

i.e. F is a smooth function on M" depending on a parameter A, and periodically-forced 
systems 

d X 

— =F{x,t,X) xMxMk^M", (1.2) 

i.e. F also depends periodically on the independent variable t. In gl and fJSJ we will 
describe (closely related) algorithms for approximating the invariant tori created at 
Neimark-Sacker bifurcation points of both (jl.ip and (|1.2D : Neimark-Sacker bifurca- 
tion for (|1.2p is defined by Assumptions 14.11 to 14.51 on pages [l2l[T6l while Neimark- 
Sacker bifurcation for ()1.1|) is defined by Assumptions l5.1l to l5.5l on pages [2314271 In f}2l 
however, we first introduce some of our ideas within the relatively simple paradigm 
case of Hopf bifurcation for (ll.ip . which is defined by Assumptions 12.11 to 12.31 on 
pages ISHU since the periodic orbits created here depend only on a single frequency. 
In contrast, the invariant tori created at a Neimark-Sacker bifurcation point have two 
independent frequencies and it is their possible resonance that creates the difficulties. 

Let (a;*, A*) e M" x M be a stationary solution of (HH]), i.e. F{x*,X*) = 0, at 
which J (a;*, A*), the Jacobian matrix of F, has n — 2 hyperbolic eigenvalues (nonzero 
real part) and a pair of purely imaginary eigenvalues. By the Implicit Function 
Theorem there is a locally unique curve of stationary solutions, parametrised by A, 
satisfying 

F{x*{X),X) = 

and the key condition for Hopf bifurcation is Assumption 12.31 on page [4l that the 
two critical eigenvalues of J (a;* (A), A) must cross the imaginary axis transversally 
at A = A*. There is then a locally unique curve of periodic orbits for (jl.ip in the 
neighbourhood of (a;*, A*). Analytical methods to investigate Hopf bifurcation are 
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contained in [51[Tni[IIl[I71[Tnj. In we show how low-order Fourier approximations of 
these periodic orbits simultaneously provide approximations to both the near-identity 
polynomial mappings which locally transform p.ip to normal form and also to the 
Lyapunov coefficient in the normal form. 

For Neimark-Sacker bifurcation of (jl.ip . we assume that u*{t) is a periodic orbit 
for A = A*, of period 2ttT*. We also assume that n — 3 of the Floquet exponents of 
u* are hyperbolic and 2 purely imaginary, the other being zero of course. Hence, by 
the Implicit Function Theorem, has a locally unique curve of periodic orbits, 

parametrised by A, and our key condition is again Assumption 15.31 on page 1241 i.e. 
that the critical pair of Floquet exponents crosses the imaginary axis transversally at 
A = A*. In contrast to Hopf bifurcation, however, we need two additional conditions 
in order to guarantee the creation of invariant tori at (u* , A*): 

• the no strong resonance Assumption 15 .41 on page 1261 so that torus bifurcation 
rather than subharmonic bifurcation is generic [14j : 

• the real Lyapunov coefficient is nonzero in (j5.26p . which is equivalent to the 
parameter A moving away from A"^ at leading order in Assumption 15.51 on 
page[27l 

Analytical methods to investigate Neimark-Sacker bifurcation are contained in |T51[ni 
134] . In H5.ll we first show how Assumption 15.41 permits the computation of low-order 
Fourier approximations for our invariant tori. The information contained in these 
low-order approximations is then used, together with Assumption 15.51 to construct 
higher-order Fourier approximations in H5.31 

Neimark-Sacker bifurcation of (|1.2p is similar. We assume that u*{t) is a periodic 
orbit at A = A* and also that n — 2 of the Floquet exponents are hyperbolic, while 2 
are purely imaginary. Hence, by the Implicit Function Theorem, (|1.2p has a locally 
unique curve of periodic orbits, parametrised by A, and our key condition is again 
Assumption 14.31 on page [131 i-e- that the critical pair of Floquet exponents crosses 
the imaginary axis transversally at A = A*. We still need the above two additional 
conditions. Assumption 14.41 on page I15l and Assumption 14.51 on page [TBI in order to 
guarantee the creation of invariant tori at (it*. A*). Analytical methods to investigate 
Neimark-Sacker bifurcation for (|1.2p are contained in [T^l E] • In SJU we again use 
Assumptions 14.41 and 14.51 to first construct low-order and subsequently higher-order 
Fourier approximations for our tori. (We have chosen this ordering for the sections 
because the absence of a zero Floquet exponent simplifies our equations, in particular 
the torus parametrisation is simpler. Hence transforming (jl.2p to ()l.ip . by adding 
time as a new state variable, is not recommended.) 

The fundamental idea behind the present paper is to use the approach in [28] , 
of which [29] is a concise version, to develop the analytical foundations of a practical 
computational algorithm for approximating the invariant tori created at Neimark- 
Sacker bifurcation points. [28j actually proves the existence of invariant tori in two 
ways: 

• constructing tori invariant with respect to the vector field, which is the ap- 
proach used in the two key books [TU [3D| on the left-hand side of Figure 11.11 

• constructing curves invariant with respect to the Poincare map, which is 
the approach used in the two key papers [HI on the right-hand side of 
Figure 11.11 

We do not wish to depend explicitly on the trajectories of ()l.ip or ()1.2p and so we 
follow the vector field approach; our algorithm being based on Fourier approximation 
and Floquet theory, in particular Floquet exponents, as introduced in [22'. (Thus to 
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Fig. 1.1. Key references for Neimark-S acker bifurcation 



appreciate fully the present paper, a familiarity with the left-hand side of Figure ll.ll is 
recommended.) Hence we emphasise that in this paper our concern is with invariant 
tori as manifolds, and we neither consider the trajectories thereon nor the stability 
of the tori. (Such questions may be answered at the post-processing stage, and are 
dealt with in several of the references.) As far as we are aware, the invariant manifold 
approach in ^28j has not been developed further for Neimark-Sacker bifurcation in the 
literature, and has certainly not been combined with the Fourier approximation ideas 
in [30| . On the other hand, there has been quite a lot of related work on the invariant 
curve approach, and we refer to [17] for details and references. In the present paper, we 
first see, in i J2.1[ how straight-forward it is to construct Fourier approximations for the 
periodic orbits created at a Hopf bifurcation point, and then attempt to generalise this 
algorithm for Neimark-Sacker bifurcation. The latter has two additional difficulties: 
coping with possible weak resonances and implementing efficiently the ideas behind 
centre manifold reduction and normal form transformation. (Sections [21 [4] and [5] have 
been deliberately written to be as similar as possible, both as an aid to the reader and 
so that the key differences stand out more clearly.) Finally, we mention that [30] is not 
explicitly concerned with Neimark-Sacker bifurcation, merely with the continuation 
of invariant tori using a Fourier-Galerkin approach. In [28) . however, it has already 
been shown how Neimark-Sacker bifurcation can be reduced to this case, and the 
constructive approach in '30] is much more relevant to us than the uniform norm 
analysis based on elliptic regularisation and smoothing operators employed in . 

2. Hopf bifurcation. We start with our two basic conditions. 
Assumption 2.1. {x*,X*) e R" x M is a stationary solution of (jl.ip . i.e. 



Fix*,X*)=0. 
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Assumption 2.2. (a;*, A*) is a Hop f bifurcation point for p.ip . i.e. 3E* e M"^" 
and non- singular P* G M"^" sttc/i i/iat 

J(a;*, A*)P* = P*E*, 

where 

E* = ~ 1 , , with ={ \ ^ CO* > 

I E* y E* £ k("-2)x("-2) / 

and E* having no eigenvalues on the imaginary axis. 

This invariant subspace decomposition is a stronger assumption than required for Hopf 
bifurcation (the standard case of E* having no eigenvalue which is an integer multiple 
of \Ld* being considered in [22]) and is chosen so that this section agrees more closely 
with S]4]and [JS] From these two assumptions, the Implicit Function Theorem tells us 
that there is a locally unique curve of stationary points, smoothly parametrised by A, 
and satisfying 

F(a;*(A),A) =0. (2.1) 

The invariant subspace decomposition may also be smoothly continued locally, and 
so we have 

J(a;*(A),A)P(A) = P(A)E(A) P : R ^ R"^", E : R M"^", (2.2) 
where P(A) is non-singular and 

/e(A) o\ E:R^R2x^ 

\0 E{X)J E:R^R("-2)x(n-2) 

with 

_ f an{X) -ar{X) \ a, ^ ^ 

a, 



^^^^ - > «.(A) a,(A) 



Finally, the key transversality condition must also hold. 

Assumption 2.3. Transversal crossing of critical eigenvalues, i.e. 

al = d«(A*) ^ 0. 



2.1. Crandall— Rabinowitz formulation. We seek periodic orbits of (jl.ip in 
the form 

x*{X) +eP{X)z{d) z:S^^W (2.3) 

with unknown z, and also with unknown frequency a; G M. e plays the role of a 
small amplitude parameter, upon which the unknowns z, X and uj depend. Thus the 
periodic orbits must satisfy 

f(x*{X) + sP{X)z{e), a) - iu^{x*{X) + eP(A)^(0)} = (2.4a) 
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and the scalar amplitude and phase conditions [TUl [21 [121 H2] 

_ f{zi9),a*ie)) 



'^^'^^[{z{eia*{e)))=^'-^ (2.4b) 
here 

a*{9) EE (cos 61, sin 61, 0,...,0)^, a*{e) = (- sin6l, cos6', 0, . . . , 0)^, ei = (1,0)^, 
with the inner-product defined by 

{wi{6),W2{9)) = ^ [ wi{9)-W2{e)de (2.5) 

for Wi,W2 ; §^ ^ M'* and s > 1. In order to apply the Implicit Function Theorem 
to (|2.4p . we must first eliminate the curve of stationary solutions: thus the Crandall- 
Rabinowitz formulation (as used in [6] for the bifurcation of non-trivial stationary 
solutions) writes 

0(2(0), A; e) = ^P{\)-^f(x*{X) + eP{\)z{e),\j 
and solves (|2.4p in the form 



= T{z{e),X,uj;e) = { ' J dO ' ' (2.6) 

7(2) - ei. 

Thus, using (|2.H) and (|2.2p . we can expand G in the form 



G(2,A;£j = E(A)2 + ^eP-iGp(2;Aj Gp : K" x K t-^ M", (2.7) 

p>2 

the n components of Gp being homogeneous polynomials of degree p in the n com- 
ponents of z with coefficients depending on A. (Here and later we display important 
functions and mappings in this way; with the understanding that the sum is limited 
by the available smoothness.) At e = 0, (|2.6p has the solution 

z{e) ^a*{9), A = A*, w = w* 

and the linearisation of (|2.6p about this solution is non-singular since a simple Fourier 
analysis (using the properties of E*) shows that 



z{e) + AE(A*)a*(6l) - ua{e) = f{e) 
7(2) =0 

implies the existence of a constant Cc > such that 

max{||2||^i,|A|,|c^|}<C£||/|U2. 

(Here we use the standard spaces/norms of periodic functions j30j . based on the inner- 
product (j2.5p .) Hence the Implicit Function Theorem, relying on a Newton-chord 
iteration for constructing solutions of (j2.6p from the starting value 

z'-°H9) = a* (9), A(") = A*, J"^ = Lu*, 
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gives the following result. 

Theorem 2.1. For all \e\ sufficiently small, (12. 6p has a locally unique solution 

X*{e), Lu*{e) and z*{e]e). 

It can he written as an expansion in powers of e \14^ , i.e. 



p>i p>i 
p>i 



(2.8) 



where only depends on the even Fourier modes 0, 2, 4, . . . , 2p and 2;2p+i onZy de- 
pends on the odd Fourier modes 1, 3, 5, . . . , 2p+l. The amplitude and phase conditions 
force 

{z*^^,ie),a*{e)) = = {z*^^,{e),a*i9)) . 

p.Sp can also be expressed in Fourier modes, i.e. 

z*{e;e) EE a*{e)+al{e) + ^ <,(£) cos m0 + 6*„(£) sin m6', (2.9) 

m> 1 

where 

a^le) has terms e,e'^,e^,... 

a\{£),b\{£) have terms e^,£^,e^,... 

m>2 a*„X'^),h*J,£) have terms 6"-^ £"+3, . . . . 

Again, the amplitude and phase conditions force 

a^(£) • ei +6*(£) ■ 62 = and a^(£) • 63 - 6^ (£) ■ ei = 0. (2.10) 



In practice we can construct accurate approximations to our periodic orbits by 
computing A, lu and a finite Fourier series 



M 



Zm{6) = o.*{S) + Qo + ttra COS mO + bm sin mO 



m—l 



which solve the Galerkin equations for (|2.6p : i.e. 

S„g(^,,((?),A;£) -u;^z,,{9) = 
■jIzm) - ei = 0, 

where 5m ■ Li^ ^ LP' is the operator which performs the Fourier series truncation. 
Thus we have the usual approximation result in terms of the decay of the Fourier 
modes in ([2J1) . 

Theorem 2.2. For all |£| sufficiently small, ()2.1ip has a locally unique solution 

^M(e): ^m(£) and 

(2 12) 

^f,(0; £) = a*(0) + a^(£) + ^ a^,(£) cosrn^ + 6^(£) sinme, 

771 — 1 
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which satisfies the error bounds 



max{|K,(.;£) - S„^*(.;£)||^i, |Af,(e) - A^(e)|, |c.f,(e) - 

<CM\{\-S,,)z*{.;e)\\H^. 

(In this paper, we will not be considering any superconvergence phenomena.) 

The Fourier approximation in Theorem 12.21 has no restriction on the size of M 
and, similarly to (|2.8p . it can also be written as an asymptotic expansion in powers 
of e. Thus instead of considering the approximation error for fixed e as M increases, 
it is also possible to consider this error for fixed small M as £ — > 0. In ^2.2\ we will 
make particular use of the approximation for M = 3 , i.e. 

t^f(e) and 

3 ('2 13) 

z^{e;e)=a*i0) + a^{e) + Y,cirni^)cosm9 + b^^ie)smm0, 

m—l 

where, as in (j2.10p . the amplitude and phase conditions force 

af (e) • ei + bf (e) • 63 and af (e) • 63 - bf (e) ■ ei = 0. (2.14) 



Corollary 2.3. Comparing (I2.13P with the exact solution in (12. 9p gives the 
errors 



m = 0,2 
rn = 1, 3 



|A*(£)-Af(£) 


= 0(e4). 






= 0{s^), 






= O(e'), 






= Oie') 


TO = 2, 




= 0(8% 






= 0{s') 


TO = 1, 3 



2.2. Normal form and its Fourier approximation. Our algorithm for Hopf 
bifurcation in tJ2. II requires neither reduction to the centre manifold nor transforma- 
tion to normal form. For Neimark-Sacker bifurcation, however, these two procedures 
have to be implemented approximately in order to cope properly with possible weak 
resonances. Thus we now choose to illustrate our later approach in the present rela- 
tively simple setting. 

Instead of carrying out the standard theoretical centre manifold reduction and 
normal form tranformation [9l[T7l[25], we adopt the operational approach in |!5l [T2l[T4l 
117] and construct the necessary transformations in order to simplify the key equation 
dm), i.e. 

G(z(.),A;e)-4.(.) = ^^^^^^ 
7(2;) - ei = 0. 

By introducing 

^ f^T^ ~t\ ^ ^.^^ z (E M."^ and z e 
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we first write 

G(^,A;e) 




G 
G 



l2 



We then aim to simplify the lower terms in G and G as much as possible by con- 
structing suitable mappings 



h 



and h 



where h is a homogeneous quadratic polynomial with A-dependent coefficients and 
h is the sum of homogeneous quadratic and cubic polynomials with A-dependent 
coefficients, to define the near-identity transformations 



z = 



y + \h{£y; A) 
y 



and 



e[[yl+ yl] ao(A) + [yl - yl] a2{\) + 2yi^2b2(A)} 

[yi [yl + yl] Si (A) + m [yl + yl] bi(A) 
+yi [yl - iyl] S3 (A) + m [^yl - yl] 63 (A)} 



y + \h{ey] A) 

y + £ { [y? + yl] So(A) + [yl - yl] S2(A) + 2yiy2b2(A)} 



(2.16a) 



(2.16b) 



The homogeneous polynomials are given the above bases in order to link up with 
the Fourier coefficients through elementary trigonometrical identities, as the table in 
Figure 12.11 shows. (Of course, by writing our Fourier series in exponential form, this 



cos 6 


yi 


cos 26* 


y'i-yl 


COS 36' 


yi{yi-iyl) 


sin6' 


y2 


sin 26* 


2yi2/2 


sin iO 


y2{3yi ~ yi) 



Fig. 2.1. Linking Fourier modes and polynomials 

correspondence is simpler; but we do not wish to give the impression that complex 
arithmetic is necessary.) Thus we see how (through y\ + yl = 1) the resonant cubic 
terms, the null-space of the adjoint of the homological operator in the usual normal 
form computations |121 125] being spanned by 



[yl + yl] 



-2/2 
2/1 



(2.17) 



and [yl + yl] 

appear through these identities, and how we must have the restrictions 

ai(A) • ei +Si(A) • 62 = and ai(A) • 62 - bi(A) • ei = (2.18) 
in the definition of h. Under these near-identity tranformations, (|2.15p becomes 

(2.19a) 

^ " (2.19b) 

(2.19c) 



G^(y(0),y(0),A;e) -c.^y(0) =0 

G^(y(e),y(^),A;e) -c.^y(0) -0 
7(y) - ei = : 
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f xM"-2 xR>->M"-^ 



the two mappings 

X X K X K ^^^^ x K"-^ x K x M IR"-^ 
capable of being expanded, like (|2.7p . in the form 

(y, y, A; e) - E(A)y + ^ eP-^Gp (y, y; a) g], 

p>2 

G^ (y, y, A; e) = E(A)y + ^ e^-iG^ (y, y; a) G^ 

where the components of G^ and Gp are homogeneous polynomials of degree p in the 
components of y and y, with coefficients depending on A. Now we choose h and h so 
that the lower terms in G and G may be simplified in the following way: 

• h forces the coefficients of the quadratic terms for y in Gj to be zero 

• h forces the coefficients of the quadratic terms for y in G2 to be zero and the 
coefficients of the cubic terms for y in g\ to take the form 

[y! + yl]B{X)y, where B(A) ^ ^^"jjj 'j^f^^^^ (2.20) 

and we call the elements of this matrix Lyapunov coefficients. 
(I.e. after transformation, only a multiple of the resonant cubic terms (|2.17p remains.) 
After this simplification, if we now insert 

y{e) = a*{e) (= (cos6l, sin^)^) and y{e) = (2.2fa) 

A = A^-s^^ and ^ ^ ^ ..Mmin^Mm^n (2.21b) 
anXX*} Q!k(A*) 

into the left-hand side of (|2.19p . we can easily see that the remainder is 0{e^) for 
(|2.19ap . O(e^) for (I2.19bp and zero for (|2.19cp . Consequently, by transforming (|2.21ap 
back through (|2.16p . i.e. 

z{e) = a*{e) + ^h{ea*{e):\) and z{9) = ^h{ea*{0); A), (2.22) 

we obtain an asymptotic solution for (|2.15p . Since Theorem 12.11 already displays such 
a solution, i.e. A*(£), w*(e) and 

this must match with (|2.21bp and 11221) ■ Thus we obtain 
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and, through ^TF^, 

z*{e;e) = a*(6l) +£|ao(A*) + a2(A'^) 008 26* + 62(A*)sm 26*1 
+ |ai(A*)cose'+Si(A*)sin6i 

+a3{X*) cos 36* + 63(A*) sm36i| + 0{e^) 
z*{e; e)=e {5o(A*) + 52(A*) 008 26* + 62(A*) sin 26*} + 0{e^). 



(2.24a) 
(2.24b) 



FinaUy, by comparing (|2.24[) and (|2.8[) . we see that the coefficients defining h{.;X*) 
and h{.; X*) in (|2.16p are given exactly by the coefficients of the Fourier modes in the 
£2(0) and zl{9) terms of z*{6; e) and in the Z2{0) term of z*{9; e) for (|2.8p . Moreover, 
by comparing (I2.23P and ()2.8p . we also see that the Lyapunov coefficients /3r, (A*) and 
Pi{X*) in (I2.20p are given exactly by 

P^{X*) ^ -a^{X*)X*2 and (3,{X*) ^ uj* ~ a,{X*)X*2. 



To calculate the expansion in (j2.8p . however, requires (through G) explicit knowl- 
edge of the second and third derivatives of F, so it is practically much more convenient 
to approximate not only the coefficients of h{.; A*) and h{.; A*) but also the Lyapunov 
coefficients Pn^X*) and /3i(A*) by using instead the M = 3 Fourier approximation in 
(I2l3l) . i.e. Af (e), cjf (e) and 



z^{6; e) = a^(e) + ^ 2^(e) cosmO + b„(e) sin rnfl 



m— 1 

3 



zf^{9; e) = a^(e) + ^ o'm(£) cos m9 + b„(e) sin rn0. 



Theorem 2.4. Using the asymptotic error results in Corollary on 
our practical approximate formulae become 

/3h(A*) = -a,{X*)^^^^l^ + 0(6^), 



page\]\ 



m = 0,2 


am{X*) = i2,f,(e) + 0{e^) 


5,„(A*) = i5^(£)+0(£2) 


m = 


6,„(A*) = i6r,(e) + 0(£2) 


bm(A*) = iC(e)+0(e') 


m = 1, 3 


S™(A*) = ^af„(£) + 0(£2) 


6™(A*)-^S^(£) + 0(£2) 



We conclude by emphasizing how the M = 3 Fourier results will be used later in 
Neimark-Sacker bifurcation. For a chosen value of e, we can easily compute z^{9; e), 
X^{e) and i^3"(e) from Theorem 12.21 the two scalar outputs then give us approxima- 
tions for the Lyapunov coefficients /3r(A*) and /3i(A*), while the Fourier components 
of z^{9; e) provide approximations for the coefficients of the polynomials h{.; A*) and 
h{.; A*). With regard to Hopf bifurcation itself, the above approximate formulae may 
be regarded as alternatives to those suggested in [101 HZl US] ■ 
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Lyapunov coefRcients 



Coefficients of h{., A*) 




Quadratic coefficients of h{., A*) Cubic coefficients of h{., A*) 

10 ' ^ 10° 




Fig. 2.2. Errors in approximation from Theorem\2.4 



2.3. Numerical results. Now we illustrate the above approximations on the 
well-known Lorenz equations [HI UHl [IZl HH] 

ii — a(x2 — x\) X2 = Axi — X2 — xiXs is = xiX2 — bx3. 

a and b are regarded as fixed parameters and A is our continuation parameter. For 
a > b + I there is a subcritical Hopf bifurcation from the stationary solution curve 

a;*(A) = (^b{\-l), v/6(A-l),A- l)^ for A > 1 at A* = ^('^ + 

V / (7 — 0—1 

6(a-& - 1) 
2(cj*2 + (cr + 6+l)2)- 



where uj* = \/b{X* + a), E* = -{a + 6+1) and d* = 



We use the standard parameter values (ct, fe) — (10, |), which gives A* ~ 24.74, and 
Figure [221 displays the error for the approximations contained in Theorem 12.41 Thus 
the O(e^) convergence is verified. 

3. Computational Floquet Theory. Floquet theory enables us to transform 
linear, periodic ode's to constant-coefficient form: this both simplifies the analysis and 
leads to much more efficient approximation by Fourier methods. A detailed discussion 
is contained in j22], here we only describe concisely the results that are required. If 
the linear, periodic system we wish to solve is 



dv 



(6) + A{0)v{9) ^ f{e) v,f -.S^ ^ W\ A : 



(3.1) 
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then our Floquet-values and Floquet-vectors solve the corresponding eigen-problem 

-P{0) + A{0)P{0) ^P{9)E E e M"^", P(6') e M"^". (3.2) 

In general, to avoid the explicit use of complex arithmetic, it is necessary to work 
with both periodic (3^+) and anti-periodic mappings, i.e. if s : R i— > M then 



s ey± iff s 

Thus, more specifically, E and P in 



E = 
P = 



E+ 
E_ 



- 2tt) = ±s{9) ye e R. 

have the form 

E+ e M"+^"+, E_ e R" 



I.e. 



P±(6') G R""""* and P±(6' + 27r) = ±P±(6') V6' G R, 

for some n+, rt_ G Z with ri+ + 7t,_ = n and < < n. Then to transform 
constant-coefficient form, we transform v, f to Floquet variables 



v{0) ^ P{e)w{9) 

m - p{e)g{e) 



P+{0)w+{0) 

p^{9)w^{e) 

p+(%+(^) 

P-{0)g_{0) 



and hence arrive at the two equations 



g+ G 3^, + , g_ G y_ 



to 



-«;_(0) + E_t«_(0)=g_((?) 



in 3^"^ and 3^"" respectively. 

Finally, we emphasise that n± are not in general unique, but can always be chosen 
so that the imaginary parts of the eigenvalues of E± {the Floquet exponents) lie in 

V 2' 2^ 



N.B. For simplicity, we shall assume in |j4]and fj5]that 
n+ = n and n_ = 0. This is briefly commented on in iJS] 



4. Neimark Sacker bifurcation for periodically-forced systems. We may 

assume that the forcing in (|1.2p is 27r-periodic, and emphasise this by using as the 
independent variable from now on, i.e. (jl.2p becomes 



^ie)^F{v,0,x) 

We start with our two basic conditions. 
Assumption 4.1. At \ ^ X* , v* : 



d 



I" X §1 X 



I" is a periodic orbit of (|4.ip . i.e 



(4.1) 



F(t;*(0),0,A*)--t;*(^^)=O 



V6I G 
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Assumption 4.2. If we apply the Floquet theory in ^to 
then (|3.2p becomes 

where E* £ R"^" and P* : ^ K"^" with P*(6') non-singular \id £ §S and we have 
the invariant subspace decomposition 

E* = ~ ^ wif/i E* = , ^ CO* > 

and E* having no eigenvalues on the imaginary axis. 

The Implicit Function Theorem then gives us a locally unique curve of periodic orbits, 
smoothly parametrised by A, and satisfying 

F{v%e;X),e,\)-^v*{0-X)^O. (4.2) 

The Floquet variables in the invariant subspace decomposition can also be smoothly 
continued locally, and so we have 

P ■ X R I > R"^" 

JK(0;A),0,A)P(0;A)-P(0;A) = P(^;A)E(A) ' ^ (4.3) 



where P(0; A) is non-singular and 



E(A)^ 

with 

E(A) 



E 



E 



p2x2 



E(A) ^0 \ 

, E(A)i E:R^R("-2)x(«-2) 



"/(A) ^^(A) ) a 



Finally, the key transversality condition must also hold. 

Assumption 4.3. Transversal crossing of critical Floquet exponents, i.e. 



4.1. Crandall— Rabinowitz formulation. To start with, we attempt to mimic 
our approach for Hopf bifurcation in i i2.1l and seek invariant tori of (j4.ip in the form 



with unknown 2:, satisfying 



2 : §1 X §1 



(4.4) 



F[v*{e-\)+eP(d;\)z{e,(l)),6,\ 

d r 



— v''{6-X)+eP{6-X)z{e,< 



del 



d 



v*[e-\) + eP{e-\)z{e, 



(4.5a) 
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for some unknown w G R. Thus we are no longer following trajectories of p.2p . but 
characterising the invariance of (|4.4p by insisting that the vector field must lie in its 
tangent space |20[ ?IT\ . (|4.4p and (I4.5ap are based on a particularly simple choice of 
parametrisation for our torus, and we shall see in ij4.4l that more subtlety is required 
later. The present choice, however, is the natural analogue of Hopf bifurcation (with 
u playing the role of frequency in 4>) and enables us to approximate the normal form 
in §4.21 Of course, we also require the scalar amplitude and phase conditions 



7(^)^ 



= ei. 



(4.5b) 



where 



a*((?!)) = (cos (/>, sin 0, 0, . . . ,0)"^, a*(0) = (~ sine/), cos 0, . . . ,0)"^ 



with the inner-product defined by 



[2n] 



2tt p2tt 



■ 'W2[ 



I de d(j). 



To attempt to apply the Implicit Function Theorem to (j4.5p . we must first elim- 
inate the curve of periodic orbits: thus the Crandall-Rabinowitz formulation writes 



g(^(0, 0), a, 6; e) = iP(^?; A)-i |f(^*(0; A) + eP(0; X)z{0, ^),e 



,A 



and solves (14.51) in the form 



O^T[z{e,c^),\u,e-e) 



G z( 



, A,0;e) - Tr^(^,' 



de 

l{z) ~ ei 

Hence, using (|4.2p and (|4.3p . we can expand G in the form 



dz 



g(^,A,0;£) =E(A)^ + ^£f-iGp(^; 



A, 6* 



Gp : M" X M X M", 



(4.6) 



(4.7) 



p>2 



the n components of Gp being homogeneous polynomials of degree p in the n compo- 
nents of z with coefficients depending on A and 9. At e = 0, (|4.6p has the solution 

2(61,0) =a*(</>), A = A*, uj^uj* 

and the linearisation about this solution is 



AE(A*)a*(</)) -tja*(0) 



(4.8) 



l{z)- 



Unlike Hopf bifurcation, however, there is no guarantee that the linearisation 
is non-singular since 







-1^ 




) 


86 



d 



= 
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may have the solution z{9, (f>) = {zi{0, 0), Z2{9, 4>))^ with 

zi (61, 0) + iZ2 {9,(1)) = e'('^''+™'^) £, m e Z. 

This occurs if lo*{1 ~ m) — £, and so in particular for {£,m) = (0,1); but this is 
the same as for Hopf bifurcation and again compensated for by the scalar unknowns 
X,uj and the scalar conditions 7. Now, however, there is a difficulty whenever u* is 
rational, i.e. the resonance situation 

(4.9) 



1 — m 

One theoretical answer to this problem is to assume that lo* is not only irrational, 
but also satisfies a Diophantine condition implying that it is badly approximated by 
rationals; i.e. {l,rn) must be large in order to approximately satisfy (|4.9[) . This 
is the approach used in KAM theory [23], but here we can make a pair of simpler 
assumptions. 

Assumption 4.4. No strong resonance, i.e. 



(Here we must remember that our form of Floquet theory in ^ enforces the bound 
< u* < ^.) This assumption is required because oj* being rational is also a necessary 
condition for subharmonic bifurcation of p.2p to occur [T3]. For rational uj* with 
denominator > 5, the torus bifurcation is generic; while for w* = ^, the subharmonic 
bifurcation is generic. (For u* — j, the relative size of certain parameters determines 
whether torus or subharmonic bifurcation occurs [14[ 133] , but for simplicity we omit 
this case.) 

Theorem 4.1. Under Assumption \4-4[ ^2 can expand ()4.6p in powers of e and 
construct an asymptotic solution 

X^{e), Lu^{e) and z'^ {9 , (j); e) 

up to and including the term, i.e. 

X^{e)^X*+e'Xl c.^(e)^c.*+eV, 

POP (4-10) 
z''{9,(j);s) = a*{(j))+ez^{9,^)+e^zi{9,(b); 

where Z2 only depends on the Fourier (j)-modes and 2 and zf only depends on the 
Fourier (j)-modes 1 and 3. The amplitude and phase conditions force 

{{zi{9, 0), a*(0))) = = {{zi{9, 0), d*(0))) . 

(|4.10p can also be expressed in terms of Fourier (j)-modes, i.e. 

3 

z^{9,(j);e) EE a*{(j)) +a^{9;e) + ^ a,^(6l; e) cos m(/) + 6^(6*; e) sinrnc/), (4.11) 



where {9; e) , a|'(0;e), ^2 (^i^) aree-terms and af {9;e), {9;e), a^{9;e), 63 {9;e) 
are s'^ -terms. Again, the amplitude and phase conditions force 

{af{9; e), e,) + {bf{9; e), 62) - = {af{9; e), 63) - {bf{9; e), ei). (4.12) 
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Assumption 14.41 is also sufficient to approximately solve (|4.6p with M — 3 Fourier 
0-modes; i.e. 



Soo,3G(z3(^, 0), A, e; e) - ^z^iO, 0) - ^^z^iO, 0) = 
7(23) - ei = 0, 



(4.13) 



where the operator S00.3 '■ ^ performs the Fourier (/)-mode truncation. 

Theorem 4.2. For \£\ sufficiently small, (|4.13|) has a locally unique solution 

iFf\ , .Fi 



A (e), (e) and 
3 

z^{9, 0;e) = a*(0) + a^(6'; e) + ^ am(^;^) cosm0 + b^j(6'; e) sini 



(4.14) 



(4.15) 



j4s m ()4.12p . t/ie amplitude and phase conditions force 

(af (0;e),ei) + (bf (0;e),e2) = - (af (0;e),e2) - (bf (0; e), ei). 
Comparing \A.\1\ and (|4.14[) . as m Corollaru \2.!A gives the errors 

|A^(£) - A^(e)| = 0(^4), \oj^{e) - oj^ {e)\ = 0(e4), 

m = 0,2 ||a,f,(0;e)-a^,(0;e)|| -0(e3), ||b^(0; e) - b^^(0; e) || = ©(e^) m = 2, 

m=l,3 \\ai{e-e)~a^,{e-e)\\^0{e^), b^{0; e) ^ b^iO; s) = 0(e*) m=l,3. 



We can now state our second condition, which may be expressed in several equiv- 
alent forms. 

Assumption 4.5. Nonzero real Lyapunov coefficient, i.e. 



Af ^i[A^]"(0)^0. 



Since A^(e) — A* and A^(£) — A* have no 0(e) term, Assumption 14.51 forces A^(£) 
and A^(e) to move away from the critical value A* for small e 7^ 0. Together with 
Assumption 14.41 it also shows that aj,(A^(e)) and aR{\^{e)) move away from zero 
for small e 7^ and therefore permits merely the no strong resonance condition in 
Assumption 14.41 (This pair of assumptions has its analogue for Hamiltonian systems 
[21].) We shall see later, in (|4.25p and Theorem l4.31 that Assumption 14. 51 is equivalent 
to a real Lyapunov coefficient being nonzero. 

4.2. Normal form and its Fourier approximation. In order to cope with 
possible weak resonances, we need to reduce our equations to an approximate normal 
form. Our algorithms in ij2] for the existence, uniqueness and Fourier approximation 
of periodic orbits created at a Hopf bifurcation point required neither reduction to the 
centre manifold nor transformation to normal form: for Neimark-Sacker bifurcation, 
however, these two procedures have to be implemented approximately and in this 
subsection we follow the strategy in ij2.2l 

Our aim is to simplify the key equation (j4.6p . i.e. 



G 



£ - 



d_ d_ 







(4.16) 



7(2) - ei = 0. 
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By again introducing 

z = {z, zf , with zeR^ and z G R"^^ 

we can write 

G{z,z,X,9;e' 
G(z,z,\e;e] 

We then construct 



G 



G : IE 
G : 



En — 2 



where his a homogeneous quadratic polynomial with {9, A)-dependent coefficients and 
h is the sum of homogeneous quadratic and cubic polynomials with {9, A)-dependent 
coefficients, to define the near-identity transformations 



2 - y + ^h{ey; 9, A) 

= y + e{[yl+yi] MO, A) + [yf - yl] MO, A) + 2M2b2{d, A)} 
+ {yi [fi + yl] ai{0, A) + h [yl + yl] MO, A) 
+yi [yl - 3yl] Md, A) + h [3yl - yi] A)} 



(4.17a) 



and 



z = y + ^h{ey;9,X) 
= y + £{[yf+yi] M^, A) + [yf - yi] M9, A) + 2yiy2b2{9, A)} . 
As in (|2.16p , we must also have the restrictions 

{MO, A), ei) + (bi{9, A), 62) = = {MO, A), ^2) - {MO, A), ?i) 



(4.17b) 



(4.18) 



in the definition of h. Under these near-identity transformations, (|4.16p becomes 



G^(y(0,0),y(0,</)),A,0;< 
G^(y{9,<l>),y{9,cP),X,9;^ 



e - 



s - 



d_ d_ 

d9^'^d(j) 
d d 

h UJ 

89 86 



y{0,4>) = Q 
y{0,^)^^ 



2 „„j _ „2 X jjn-2 X M X §1 X 



7(y) - ei = : 

the two mappings 

: R2 X ]R"-2 X M X §1 X R and G 

capable of being expanded, like (|4.7p . in the form 
G^(y,y,A,0;e) = E(A)y + ^ ef-iG],(y, y; A, 0) gJ, 

p>2 

G^(y,ll,A,0;e) =E(A)y-F5]£f-iGp(y,y;A,0) gJ, 



(4.19a) 

(4.19b) 
(4.19c) 



p>2 
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where the components of and Gp are homogeneous polynomials of degree p in the 
components of y and y, with coefficients depending on A and 9. Now we choose h 
and h so that the lower terms in G and G may be simplified in the following way: 

• h forces the coefficients of the quadratic terms for y in G2 to be zero; 

• h forces the coefficients of the quadratic terms for y in G2 to be zero and the 
coefficients of the cubic terms for y in G3 to take the form 

[y! + yl]B{X)y, where B(A) ^ ^^"^j^j 'j^f^^^^ (4.20) 

and we again call the elements of this matrix Lyapunov coefficients. 
(I.e. after transformation, only a multiple of the resonant cubic terms (|2.17p remains.) 
After this simplification, and under Assumption 14.41 if we now insert 

y{e,(t)) = a*{(j)) (= (coacj), sin cjjf) and y(6l,0) = (4.21a) 

X y* -2 /^h(A-) , , ^2 «h(A-)/3.(A-)-A,(A*)/3h(A-) 

A = A — e - — 7—^ and uj — oj +e : — 7—^ (4.21b) 

into the left-hand side of (|4.19p . we can easily see that the remainder is 0{e^) for 
(|4.19ap . O(e^) for (|4.19bp and zero for (|4.19cp . Consequently, by transforming (|4.21ap 
back through (|4.17p . i.e. 

z{0,(l>) ^a*{<j)) + ^h{ea*{(P);e,X) and z{d,(j)) ^ ^h{en*{(j));d,X), (4.22) 

we obtain an asymptotic solution for (|4.16|) . Since Theorem 14.11 already displays such 
a solution, i.e. A^(£), uj'^{e) and 

z (6^,0, £) = _E 

\z {e,(f>;e) 

this must match with (|4.21bp and (|4.22p . Thus we obtain 
A^(£)-A'^-e2^"^^*^ 



aniX* 



E,,_ 2«H(A'^)/3i(A*)-a,(A*)/3,(A*) 



(4.23) 



d«(A*) 
and, through (|422l), 

£^(6*, (/);£) = a*(0) +e|ao(6l. A*) + 02(6', A*)cos2(/)+S2(6',A*)sin20| 
+ |ai(6l,A*)cos0 + bi(6',A*)sin0 

+03(6*, A*) cos 30 + 63(61, A*) sin 3(?!)| (4.24a) 
5^(6*, (/);£) = £100(6*, A*) + 02(6*, A*) cos 20 + 62(6', A*) sin 2(^1 +0(£2). (4.24b) 

Finally, by comparing ()4.24p and ()4.10p . we see that the coefficients of h{.;d,X*) 
and h{.;d,X*) are given exactly by the coefficients of the Fourier 0-modes in the 
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z^iO, (p) and 



terms of z (6, e) and the Z2 {0, 0) term in z [9, 0; e) for 



(j4.10p . Moreover, by comparing (|4.23p and (I4.10p . the Lyapunov coefficients /3r(A*) 
and I3i{\*) in (|4.20p are given exactly by 



/?h(A*) = -d„(A*)Af and (3,{\*) = ~ a:{\*)\^ , 



(4.25) 



and now we see that Assumption 14.51 is equivalent to /3r(A*) / 0. 

To calculate the expansion in (j4.10p . however, requires (through G) explicit knowl- 
edge of the second and third derivatives of i^, so it is practically much more convenient 
to approximate not only the coefficients of h{.] 9, A*) and h{.; 9, A*) but also the Lya- 
punov coefficients /3h(A*) and Pi{X*) by using instead the M = 3 Fourier 0-mode 
approximation in (I4.14p . i.e. A^(£:), uj^{s) and 



z^{6,(j);e) =aQ{6,e)+ ^ 2^ (^'i s) cos to0 + b^{6,e) sin m0 



m—1 
3 



e) cos 7710 + b^{9,£) sin mq 



Theorem 4.3. Using the asymptotic error results in Theorem \4-^ on page\ 
our practical approximate formulae become 



/3r(A*) 



Pi(A ) = 2 "H'^ ) 2 ^ '-'i^ )' 



m = 0,2 


a,„(0,A*) = iSr,(0,e) + O(e2) 


a,n{9,X*)^'-ai{e,e)+0{e') 


TO = 2 


S™(0,A*) = iC(^?,£)+O(£2) 


h^i9,X*)^^KiO,e) + 0{e^) 


TO = 1, 3 


a,^{9,X*)^^aiie,e)+0{e^) 


b^{9,X*)^^K,{9,e)+0{e^) 



We conclude by emphasizing how the Af = 3 Fourier </)-mode approximation plays 
the same practical role for Neimark-Sacker bifurcation as that described in the final 
paragraph of i J2.2l 

4.3. Numerical results. As a numerical example, we use the forced van der 
Pol equation |51 [1^1 dS] , which may be written in the form (|1.2p as 

ii = X2 + ox\ (1 — a;J/3) , ±2 = —x\ -I- A cos vt. (4.26) 

Here cr > and < < 1 are regarded as fixed parameters and A, as usual, is our 
continuation parameter: in the form (j4.ip . (|4.26p becomes 

= - {w2 + cr?;i (1 - v\lZ) } , W2 = - {-wi + A cos 9} . (4.27) 
For (T = 0, it is interesting that (|4.27p has the periodic orbit and Floquet variables 
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V 


A* 






(A*) 


«.(A*) 


/3r(A*) 


/3i(A*) 


0.86 


0.4349 





1092 


-3 


3492 


-1.9037 


-0 


1154 


-0 


1607 


0.85 


0.4536 





1271 


-2 


7032 


-1.3431 


-0 


0906 


-0 


1177 


0.84 


0.4726 





1444 


-2 


2858 


-1.0166 


-0 


0759 


-0 


0933 


0.83 


0.4917 





1614 


-1 


9948 


-0.8072 


-0 


0663 


-0 


0779 


0.82 


0.5110 





1784 


-1 


7808 


-0.6634 


-0 


0597 


-0 


0673 


0.81 


0.5304 





1953 


-1 


6173 


-0.5598 


-0 


0549 


-0 


0597 


0.80 


0.5498 





2124 


-1 


4886 


-0.4822 


-0 


0513 


-0 


0540 


0.79 


0.5692 





2296 


-1 


3849 


-0.4223 


-0 


0486 


-0 


0496 


0.78 


0.5886 





2471 


-1 


2999 


-0.3749 


-0 


0465 


-0 


0460 


0.77 


0.6079 





2648 


-1 


2290 


-0.3367 


-0 


0449 


-0 


0432 


0.76 


0.6271 





2829 


-1 


1693 


-0.3054 


-0 


0436 


-0 


0408 


0.75 


0.6462 





3013 


-1 


1185 


-0.2793 


-0 


0425 


-0 


0389 


0.74 


0.6651 





3200 


-1 


0749 


-0.2573 


-0 


0417 


-0 


0372 


0.73 


0.6840 





3392 


-1 


0372 


-0.2386 


-0 


0411 


-0 


0358 



Fig. 4.1. Neimark-S acker bifurcation points for the forced van der Pol equation 



which is useful as a starting value for continuation. (Note that the eigenvahics of E 
in (|4.28p are purely imaginary; and in fact there is a "degenerate" Neimark-Sacker 
bifurcation here, with respect to the parameter a, for which the invariant tori formu- 
lae, all at cr = 0, may be written down exactly. This is of no interest to us.) Having 
computed a periodic orbit at the value of a we are interested in, we can then fix a and 
continue in A, looking for Neimark-Sacker bifurcation points. We use the techniques 
described in [22] and, because of the form of the forcing, the periodic orbits have the 
symmetry 

v{e + Ti) = -v{e) v^eR; 

which has the important practical simplification that v{6) need only be approximated 
by odd Fourier modes. This symmetry is inherited by the Floquet decomposition in 
SJSl so that (if we use the strategy in ^22j to limit the size of the imaginary part of the 
Floquet exponents) either 

p{e + 'K) = ±p{e) and w{e + -k) = ^w{e) ve'eR. 

In Figure 14.11 we display A* for Neimark-Sacker bifurcation points at different v val- 
ues but with (7 = 4, and this may be compared with Figure 13 in [55]. (A simple 
secant iteration was used to locate the periodic orbits with purely imaginary Floquet 
exponents, so we are not using a sophisticated method to detect Neimark-Sacker bi- 
furcation points.) We want to show how some of the important scalars associated 
with the bifurcation vary with v in this example; and so we display the w*, di,(A*), 
dj(A*), /3r(A*) and /3i(A*) values at these bifurcation points, the latter pair being 
approximated as in Theorem 14.31 with e = 0.005. (Note that we have jumped across 
two points of strong resonance, where u* = \ and ^.) For these calculations we 
used M = 24 Fourier 0-modes, which reduced the size of the Fourier coefficients to 
« 10-14. 

4.4. Higher-order Fourier approximation of tori. In order to compute 
higher-order approximations for our invariant tori, we must employ a more suitable 
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parametrisation than (j4.4p . Thus we use the normal bundle of the approximate torus 

v*{9-X) +eP{e;X)a*{(j)) (4.29) 

and, in (liT^ . 

• replace y{6, (/>) with [1 + p{d, cf))] a*{(t)) for unknown p : S-'^ x S-"^ ^ M, 

• allow w : X §^ I— > M to be an unknown function . 

This links up with the invariance condition used in [21 for continuation of tori, and 
corresponds to using polar co-ordinates in the critical 2-dimensional subspace. In 
(|4.19c|) there is now no need for a scalar phase condition, and the scalar amplitude 
equation simplifies to a zero-mean condition for p, i.e. 



((p(0,0),l))=O. 
Thus our equations for p and uj in (|4.19ap decouple to become 



a*((/.) • G' ( [1 + Pie, 0)] a*{cj,),yie, 0), A, 9- e 

d 



do 



p(0,0) = O 



(4.30) 



(4.31a) 



and 



^—1— ■ ( [1 + Pie, 0)] a*icp),y{e, 0), A, 9; e) - ^{9, 



p{9,<p) 

while the hyperbolic equations in (j4.19bp remain 
G^ ( [1 + p(0, 4>)] a*i4>), viO, 4>),\ 0-e)- ^ 
The crucial leading terms in (|4.31[) are 



= 0, (4.31b) 



(4.32) 



and 



a*(0) • G' ( [1 + p(0, <P)] a*{cj,), yi9, </)), A, 9; e) 

= [1 + pi9, cl>)] aniX) + e^PuiX) [1 + p{0, + O(e') 

— ^— -^(0) • G^ ( [1 + p(0, 0)] S*(0), y(0, 0), A, 9; e) 



(4.33a) 



(4.33b) 



= a,(A) + e^f3,iX) [1 + pi9, + Oie^). 



Consequently, if we use (j4.31bp to define uj{9, (f) in terms of A, p{9, 4>) and y(0, 0) 
for the rest of this subsection, wc finally have to prove that the system of equations 



a*(0).G^([l + p(0,</))]a*(0),y(0,0),A, 

d_ 

m ^ 

G^([l + p(0,</))] a*(^),y(0,0),A, 

d 



p(0,0)-O 



{{pi9,4,),l))=0 



(4.34a) 

(4.34b) 
(4.34c) 
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has a locally unique solution {\,p,y) for \e\ sufficiently small. This is achieved in 
PSI through the iteration 



E* - 



961 



9 



(fc+i) 



90 



y ' ''{0,4') = r 



(4.35a) 

(4.35b) 
(4.35c) 



where 5X^''+^'> = A^'^+i) - A^^) and 

.W^fl,0) EE -_^i^r(0) . g\ [l + pW(f?,0)] 2-^(0), y('=)(0,0),AW,0;e) 



pW(0,0) 
f('^-)(0,0)^2e2/?^(A*)p«(0,0) 

- S*(0) • gV [l + pt'^) (0, 0)1 a*(0), y« ((?, 0), A^'^) , ( 



l + pW(0,0) a*(0),y('='(0,0),AW 



with starting values 



A(°)=A*, p(o)(0,0)=O, y 



(0), 



The key idea behind showing that these iterates remain bounded and then converge 
is to integrate (|4.35ap against ^'^'^+^•'(0,0), so that the l.h.s. becomes 



fc+i) 



2e'/3H(A*) 



d 



2£2/3,(A*) + |^^W( 



(4.36) 



Since (I4.33bp shows that the leading non-constant term in u}'^^\9,(t)) is O(e^), As- 
sumption [43] ensures that (|4.36p is a definite quadratic term for |e| sufficiently small, 
and this is sufficient for [28] to prove the following result. 

Theorem 4.4. Suppose F in (|4.ip has r > 5 continuous derivatives for (A,i;) 
in a neighbourhood of (X*,v*{0)). Then 3er > such that for \e\ < Sr (|4.34p has a 
locally unique solution 

A*(£), p%d,4>-e), y*{OA;e) 

with p*{.,.;e), y*(.,.;e) having (r — 1) Lipchitz continuous derivatives. This means 
that both y and z, z ( through (|4.17p ) have this degree of smoothness, and so, through 
(|4.4p . do the invariant tori as manifolds. 

The subtlety of Theorem 14.41 is that, in general, ~^ as r — oo; in particular, one 
cannot expect the tori to be analytic when F is analytic. 

In practice we seek an approximate solution of (|4.34p in the form 



l,rn 



\£\ < L, \m\ < M 
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with p_i _m, U-i.-m the conjugates oi pi^ra, Uim- These functions must satisfy 



a*(0) • G' ( [1 + p,,,,{e, <P)] a*{c^), y, ,,{e, 0), A, 0; e 



= 



g7 [1 + p,.,,{e, 4>)] 2*(0), ,,(0, 0), A, 0; e 



d 



d 



de 



for 1^1 < L, |m| < M and Vx e R"-^ with 



i(^e+m0) 



(4.37a) 

(4.37b) 
(4.37c) 



l + Pl,,A/(^, 0) ^ ^ 



As shown in [30], the iteration analogous to (|4.35p also converges here and gives the 
following result. 

Theorem 4.5. Under the conditions of Theorem \4.4\ (|4.37p has a locally unique 
solution 



^f.M(e), PL.M{0,(t>;£), 2/1,. m(^, 



satisfying 



\Xl,,{e)'X*{s)\ 



max^ IIp^ j,,,(., .;£) - Si,.MP*(-, •;e)||L2 > < Cmax 
l|yr,A/(-,-;£) -SL.My*(.,.;£)||L, 

We comment on the implementation of this algorithm in JJS] 



I - Sz,.m)p*(., .;e)||ffi 
' -Si.M)y*(-,-;£)||ffi 



5. Neimark Sacker bifurcation for autonomous systems. We start with 
our two basic conditions. 

Assumption 5.1. At X = A*, (|l.ip has a periodic orbit u*(t) of period 2ttT* , 
and so, under the change-of-variable 



T* and v* satisfy 



v*{e) = u*{9T*) V* : §1 



T*F{v*{e),\*) v*{e) = 0. 

d9 



Assumption 5.2. If we apply the Floquet theory in ^to 

A{e)=T*J{v*{e),X*) 

then p.2p becomes 

T''J{v*{0),X*)P*{0) - P*(e') P*{e)E*, 
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where E* G M"""" and P* : S.'^ ^ M"^" ivith P*{e) non- singular G S\ and we have 
the invariant subspace decomposition 



E* 
E* = ( E* 




1,2x2 



E* e 

E* £ ]]j("-3)x(n-3) 



with E* 



-cj* 




> 



and E* having no eigenvalues on the imaginary axis. 

The Implicit Function Theorem then gives us a locally unique curve of periodic orbits, 
smoothly parametrised by A, and satisfying 



T*iX)Fiv*ie;X),\)~—v*ie;X)=0. 



(5.1) 



The Floquet variables in the invariant subspace decomposition can be smoothly con- 
tinued locally, and so we have 



T*{X)j{v*{e; A), X)P{e; A) - P(6'; A) = P{0; A)E(A) 



where P{9; X) is non-singular and 

/E(A) 0^ 

E(A) = E(A) 

V Oy 



E 
E 



P : §1 X M K"^" 
E : M K"^", 



1^2x2 

B(n-3)x(n-3) 



(5.2) 



with 



E(A) 



_ f a^iX) -aj(A) 
aj(A) a„(A) 



an 
a, 



Finally, the key transversality condition must also hold. 

Assumption 5.3. Transversal crossing of critical Floquet exponents, i.e. 



al = d«(A*) ^ 0. 



5.1. Crandall— Rabinoviritz formulation. To start with, we attempt to mimic 
our approach in §4.11 and seek invariant tori of (|1.1|1 in the form 



v*{e-X)+eP{e]X)z{9,^ 
with unknown satisfying 



2 : §1 X 



T*{X)F[v*{e; X) + eP{e- X)z{e, q)),X 

d 



[i + e,i]-^^[v*{e-x) + eP{e-x)z{e,, 



v*{e]X) + eP{e-x)z(e, 



(5.3) 



(5.4a) 



for unknown w,?] G R. As in i i4.11 we are expressing the invariance of (|5.3p by 
insisting that the vector field lie in its tangent space; the only difference being the extra 
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unknown rj now compensating for the zero Ffoquet exponent. The parametrisation 
of the torus in (|5.4p is again special, however, since the coefficients ([1 + er/] and w) 
must be constant: as in i i4.f 1 we will need to generalise this parametrisation later in 
''^^^ Of course, we also require the scalar amplitude and phase conditions 



7i(^) ^ 



= ei. 



(5.4b) 



and the scalar phase condition 

72(^) = ((2(^^,0),e„))=O. (5.4c) 

To attempt to apply the Implicit Function Theorem to (j5.4p . we must first elim- 
inate the curve of periodic orbits (|5.ip : thus the Crandall-Rabinowitz formulation 
writes 



g(^(0, <^), a, 77, 9; e) = iP(0; Xy^ i^T*{X)F (v*{9; A) + eP{9; X)z{e, 0), A 

- T*iX)F(v*{9; A), a) - e [1 + erj] P{9; X)z 

and solves (15.41) in the form 



G{z{9,(l)),X,f],9;e) - r/e, 
d 



= T{z{9,(l)),X,r],u,9;e) = < 



(5.5) 



Hence, using (|5.ip and (|5.2p . we can expand G in the form 

g(2, a, 77, 0; e) = E(A)^ + er7Gi (2; A, + ^ e^^^Gp (2; A, 0) , 



(5.6) 



where Gp : K" x M x §^ ^ R" and the n components of Gp are homogeneous 
polynomials of degree p in the n components of z with coefficients depending on A 
and 9. At e = 0, (|5.5p becomes 



- 77e„ = 



72(2) =0 

with solution 

z{9, (f) = a*(</>), A = A*, w = uj\ 77 = 0; 

and the linearisation about this solution is 
d 



d 

E* - I— ~uj*\ — 
89 d<j) 



z{9, (t>) + AE(A*)a*(</.) - uja*{(f) - rje, 

7i(^) 
72(2)- 



(5.7) 



26 



GERALD MOORE 



Just as in ij4.11 there is no guarantee that the hnearisation (j5.7p is non-singular, 
and singularity occurs if 



(m - 1)lu* + £^0 or muj* + £ = 



for some £, m e Z. 



(5.8) 



The first occurs for {£, m) — (0, 1), but this is compensated for by the scalar unknowns 
X,uj and the scalar conditions Ji, and the second occurs for {£,m) = (0,0), but this 
is compensated for by the scalar unknown r/ and the scalar condition 72. If to* is 
rational, however, (|5.8p will be satisfied by larger integer values of {£, to); even if oj* is 
irrational, they will be satisfied "arbitrarily" closely. Thus we must impose the same 
condition as in in ^AA\ 

Assumption 5.4. No strong resonance, i.e. 



It then follows that an asymptotic solution for (|5.5p can be constructed. 

Theorem 5.1. Under As.sumption \5.4[ we can expand (|5.5p in powers of s and 
construct an asymptotic solution 

A^(£), w-^(£), r;^(e) and 2-^(6*, e) 

up to and including the terms, i.e. 



9,0;e) = a*{(j))+ez. 



(5.9) 



where Z2 only depends on the Fourier (p-modes and 2 and only depends on the 
Fourier (j)-modes 1 and 3. The amplitude and phase conditions force 

{{zi{9,cf>),a*{q^)))=Q = {{zi{9,cf>),a*{cl,))) and ((2;f (0, 0), e„)) = 0. 

(|5.9p can also be expressed in terms of Fourier (jy-modes, i.e. 

3 

z'^{e,(j)]£) = a*{(j)) +a^{e;e) + ^ a^(0; e) cos m<^ + 6^(6*; e) sin TO(/), (5.10) 

m=l 

where a^{e;e), a^{9;e), (6';e) aree-terms andaf{9;e), bf{9;e), af{9;e), bf{9;e) 
are e'^ -terms. Again, the amplitude and phase conditions force 



(af(0;e),ei) + (6f(0;£),e2) =0 
{af{9;e),e2)-(bf{9;e),e,)^0 



id (a^(0;e),e„) =0. 



(5.11) 



Assumption 15.41 is also sufficient to approximately solve (|5.5p with AI = 3 Fourier 
(/)-modes; i.e. 



S^,3G(^Z3{9,(j)),X,r],0;e^ - Tye^ 



Z3i9,(f>)=0 



71(2:3) - ei = 
72(2:3) = 0, 



(5.12) 
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where the operator Sqo.s '■ LP' ^ performs the Fourier 0-mode truncatfon. 

Theorem 5.2. For \e\ sufficiently small, (j5.12p has a locally unique solution 



A (e), (e), 77 (e) and 
3 

2^(61, 0;e) = a*(0) + a^(6';e) + ^ a,fj(6l; e) cos + 6^^(e'; e) sin 1 



(5.13) 



^(i (a^(6';e),e„) = 0. 



As in (15. lip , t/ie amplitude and phase conditions force 

(af(^;£),ei) + (5f(^;£),e2) =0 

/ \ ' 

(af(0;£),e2)-(5f(^;£),ei)=O 

Comparing (|5.10p and (|5.13p gives f/ie errors |?7^(e) — 77^(e)| = ©(ff"^), 

|A^(e) - A^(e)| = 0(£^), L^(e) - cc;^(£)| - 0(£^), 



m = 0,2 af: 



;e)|| =0(e3), 



m=l,3 \\agi9;e)-a^,ie;e)\\^Oie^), 



(5.14) 



m = 2, 
TO = 1, 3. 



We can now state our second condition. 

Assumption 5.5. Nonzero real Lyapunov coefficient, i.e. 



Af^i[Al (0)7^0 



As in ii4.1i this means that A^(e) and A''^(e) move away from the critical value A* for 
small e 7^ 0. (I5.26P and Theorem 15.31 show that this is equivalent to a real Lyapunov 
coefficient being nonzero. 

5.2. Normal form and its Fourier approximation. We follow the strategy 
in ij4.21 and construct the necessary transformations in order to simplify the key 
equation (|5.5p . i.e. 



G{z 



= 



7i(2) - ei = 
72 (^) - 0. 



(5.15) 



By introducing 



we can write 



G{z,X,'n,9;e 



z = {z,z,z) , with (z, 2, 1) G 



Glz,z,z,X,ri,6;e] 
G(z,z,z,X,ri,9]e] 
Giz,z,z,X,ri,6;e) 



G : X R""'^ xRxRxRxS^x 



X R"-'^ 

G : R2 X R""^ xRxRxRxS^x 
G : R^ X R""3 xRxRxRxS^x 
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We then construct 



h 



X §1 X 



pn— 3 



and h 



X §1 X 



where ^ is a homogeneous quadratic polynomial with (6, A)-dependent coefficients and 
h, h are the same as in H4.2[ to define the near-identity transformations 



z ~ y + ^h{ey; 6, A) with expansion as in (|4.17p 
z = y + i/j.(ey; 9, A) with expansion as in (|4.17p 
z = y+ jh{ey;0,X) 
= y + e{[yl+yi] ao{0, A) + [yf - y|] (12(6, A) + 2^1^262(0, A)} 

Now we must have the restrictions 



A), ei) + i^biie, X),e2) - - {ai{e, A), ^2) - A), ?i 

in the definition of h, and the restriction 

{ao{9,X),l) ^0 



(5.16a) 
(5.16b) 

(5.16c) 
(5.17) 
(5.18) 



in the definition of h. Finally, to compensate for (I5.18|) . it is also necessary to include 
the near-identity transformation 



77 = C + e'«(A). 
Under all these transformations, (|5.15p becomes 

G^{y{0,<P),y{0,ct>),y{d,^),KC,O;e) 

d 



(5.19) 



d 



yie,^) = o 



G'(y{e,^),y{e,^),y{0,^),XX,O;e 



{l+e^K{X)+eC)^ +UJ 



d 



80 d(j)_ 

Gt 0), 0), A, 0', e) - [C + ^^(A)] 



(l + sMA)+eC)|+^^ 

7i(y) - ei = 
l2{y) = : 



y(0,0)=O 



(5.20a) 



(5.20b) 



(5.20c) 

(5.20d) 
(5.20e) 



the three mappings 



: X M""'^ xKxKxMxS^x 

: X M""''' xKxKxMxS^x 
: m2 ^ M""^ xMxMxMxS^x 



pn— 3 
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capable of being expanded, like (|5.6p . in the form 

(y, y, y, X, (, 9; e) = E(A)y + ^ sP-'gI [y, y, y; A, 0) + C J] s^kI (y, y, y; A, 

p>2 p>l 

(y, y, y, A, C, 0- e) = E(A)y + ^ e^'^G], {y, y, y; A, fl) + C ^ s^kI (y, y, y; A, 0) 

p>2 p>l 

Gt (y, y, ^, A, c, ^; e) = E ^"^'^P (y' y' m\d)+cY. ^"K {y^ ^) ; 

where the components of both Gp, Gj, and /iTp, K^^ Kj^ are homogeneous polyno- 
mials of degree p in the components of y, y and y, with coefficients depending on A 

and 0. Now we choose ft., ft. and h so that the lower terms in G , G and & may be 
simplified in the following way. 

• ft forces the coefficients of the quadratic terms for y in g\ to be zero; 

• h, and k(A) in (|5.19p . force the coefficients of the quadratic terms for y in gJ 
to take the form 

• ft forces the coefficients of the quadratic terms for y in g\ to be zero and the 
coefficients of the cubic terms for y in g\ to take the form 

[yl+yl]%{X)y, where B(A) ^ ( ^^jjj '^pf\)) ^^''^^^ 

and we again call the elements of this matrix Lyapunov coefficients. 
(I.e. after transformation, only a multiple of the resonant cubic terms (I2.17P remains.) 

After this simplification, and under Assumption 15. 4[ if we now insert = 
together with 

y(0,</)) = a* ((/.), y(0,0)=O, y(0,^)=O, (5.22a) 

, X* ,2^^ «H , , ^2 «H(A-)/3.(A-)-d,(A-)/3^(A*) 

A = A — e — ——^ and u) — uj +e 7—^; (5.22b) 

aaiX*) Q;„(A*) 

into the left-hand side of (|5.20p . we can easily see that the remainder is 0{e^) for 
dESOil, 0(£2) for ((5:20bl) and dODcl), and zero for ffblOd^ and (jOOe)) . Consequently, 
by transforming (|5.22ap back through (|5.16p . i.e. 

z{e, 0) = a*(0) + ift(£a*(0); 0, A) 

5(0,0) = ift(£a*(^);0,A) (5.23) 
z(0,0) = U(£a*(^);0,A), 

we obtain an asymptotic solution for (|5.15p . Since Theorem 15.11 already displays such 
a solution, i.e. A^(e), Lo^{e), r}^{£) and 

/£^(^,</>;£)\ 
z^(0,0;e)^ 5^(^,<^;£) , 

Vl^(^,<^;£)/ 
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this must match with (|5.22bp and (|5.23p . Thus we obtain 



A^(e) = A* 



! /3h(A-) 
aa{X*) ' 



E _ , 2d«(A*)A(A*) - «.(A*)/3h(A*) 

U) \S) — to -f- fc 



(5.24) 



and, through ()5.23p . 



a*(0) + e |ao(6', A*) + a2(6l, A*) cos20 + 62(6', A*) sin 2^1 
+ |ai(e', A*)cos(/) + bi(6', A*)sin(/) 

+03(6*, A*) cos 30 + baCfi*, A*) sin 3<?!)| 
e |ao(6', A*) + a2(0, A*) cos 20 + 62(6', A*) sin 2^1 + 0(e2) 
e |ao(6', A*) + 02(61, A*) cos20 + &2(6', A*) sin20| + Cl(e2). 



(5.25a) 
(5.25b) 
(5.25c) 



Finally, by comparing (|5.25p and (|5.9p . we see that the coefhcients of h{.;9,X*), 
h{.;9,X*) and h{.;d,X*) are given exactly by the coefficients of the Fourier (/)-modes 



in the 22(0, (f) and z^{9, (f>) terms of {9, 0; e), the Z2 {9, (f) term in z^{9, (f>; e) and 
the i|'(^^, 0) term in z^{9,(j);£) for (|5.9p . Moreover, k(A*) = Tyf" and, by comparing 
(|5.24p and (|5.9p . the Lyapunov coefficients /3r(A*) and /3i(A*) in (|5.2ip are given 
exactly by 

/?H.(A*) = -d;,(A*)Af and P,{X*) ^ uj^ ~ d,(A*)Af . (5.26) 

Thus we see that Assumption 15 . 5 1 is equivalent to /3a (A*) 7^ 0. 

To calculate the expansion in (j5.9p . however, requires (through G) explicit knowl- 
edge of the second and third derivatives of i^, so it is practically much more convenient 
to approximate not only the coefficients of h{.; 9, A*), h{.; 9, A*) and h{.; 9, A*) but also 
the Lyapunov coefficients /3r(A*) and f3i{X*) by using instead the Af = 3 Fourier (f>- 
approximation in (|5.f3p . i.e. A^(e), uj^{s), "q^ie) and 

3 



^^(0,0;e) 



i, £) smmc? 



z^{9,(f>;e) = aQ{9,£) + ^ a^(0, e) cos to0 + b„ 

m— 1 

3 ^ 

z'^ {9 , (f); e) = a^{9,e) + ^ 5^(0, e) cos to0 + 6„ 

m— 1 
3 

£-'^(0, 0; ff) = aQ{9;e) + ^ 0^(6'; e) cosm0 + b^^{9;e) sinrnq 



7, ej smmq 



Theorem 5.3. Using the asymptotic error results in Theorem 
our practical approximate formulae are 



on page 



/3a (A*) 

MX*) 
n{X*) 



■ /\*NA(e)-A ^/ 2n 
~aa{X ) ^ h (^(e ) 

Lo^{e)-Lo* . A^(£)-A* , 
3 "AA ) 3 + C'(e ) 



77^ (e) 



+ 0(e2), 
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together with 



m = 0, 2 


m = 2 




bU0,y) = \K{0.e)+O{e^) 











and 



TO = 1, 3 



b,n{d,X*) 



^b„^{9,e)+0{e^) 



We conclude by remarking that the final comment in §4.21 applies here as well. 

5.3. Higher-order Fourier approximation of tori. In order to compute 
higher-order approximations for our invariant tori, we must employ a more suitable 
parametrisation than (j5.3p and the presence of the zero Floquet exponent in E* means 
that this parametrisation is different from (|4.29p . Thus we use the normal bundle of 
the approximate torus 



v*{9;X) +eP{9;X)a*{(j)) 



(5.27) 



and, in ([QO]) . 

• replace y{9, (/>) with [1 + p{9, cf))] a*{(t)) for unknown p : x M", 

• replace y{9, 0) by 0, 

• allow both (jj : X I— » R and ^ : §^ x i— > M to be unknown functions . 
This links up with the invariance condition used in [21] for continuation of tori, and 
corresponds to using polar co-ordinates in the critical 2-dimensional subspace. In 
(|5.20dp and (|5.20ep . there is now no need for scalar phase conditions, and the scalar 
amplitude equation simplifies to a zero-mean condition for p, i.e. 



=0. 

Thus our equations for p and uj in (|5.20ap decouple to become 



(5.28) 



a*((/)) . g' [1 + pi9, c^)] a*(^), yi9, 0), 0, A, ((0, 0; e 



d 



d 



(1 + e'K{X) + eao. 0))^ + oj{e. 



p{9,^)={) 



(5.29a) 



and 



a (0) • g' [1 + p{9, 0)] 2*(0), y{9, 0), 0, A,C(^, 0), 6- e 



l + p(0,^) 



while the hyperbolic equations in (j5.20bp remain 



G 



^ ( [1 + p{9, c^)]a*{cj,),y{9, 0), 0, A, C(e, 4>), 9; i 



{l+e'n{X)+ea9,cb))—+u;i 



u;{9,(b)^0, 



(5.29b) 



(5.30) 



y{ 
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and (|5.20cp becomes 

gU[i + p{e, cb)] a*(0), y{e, o, a, c(0, < 



[((6*, </.) + £k(A)] = 0. (5.31) 



The crucial leading terms in (|5.29p are 



a*(^) • g7 [1 + Pie, ^)] a*(</)), 0), 0, A, C(^, 0), 0; s 



[1 + p(0, 0)] a,(A) + £2/3, (A) [1 + p{d, + 0(e3) 



(5.32a) 



and 



^—l^d^ic^) • ( [1 + p(0, 0)] a*(0), y(0, <^), 0, A, C(0, 0), 0; e) 
= a,(A) + e^/3,i\) [1 + p(0, <^)]' + ©(e^), 

while in (j5.3ip we have 

Gt ( [1 + p{9, cP)] U*{4>),y{e. 0),O, A, </-), ^; e) 

= eK(A) [l + p(0,0)]' + O(e2). 



(5.32b) 



(5.32c) 



Although we needed to introduce C, through (j5.19|) in order to obtain the correct 
normal form in §5.2) it is now simpler to describe our final system of equations in 
terms of 



ri{e,^) = aO.<ty)+eK{\). 
Thus we can re-write ()5.3ip as 

G* (p(0, 0), y(0, <^), A, 77(0, 4>), 0; e) - v{0, 4>) 



= 



(5.33) 



and use (|5.33p to define j-j{6,(j)) in terms of A, p{0,(j)) and y{d,(j)) for \e\ sufficiently 
small. (Since & depends linearly on ( in (|5.3ip . and thus & depends linearly on rj 
in (j5.33p . this is particularly simple.) Similarly, we can re- write (j5.29bp as 



1 



l + p(0,^) 



a (0)-G (p{d,cl>),y{9,^),X,0;e) -uj{0,cf,)=Q 



(5.34) 



by inserting r]{9,(f>) from (|5.33p into G ; hence (|5.34p defines uj{9,(j}) in terms of A, 
p{9,(j)) and y{9,(j)). Finally, we can re-write (|5.29ap and (|5.30p as 



a*(0).G*(p(0,0),y(0,0),A,0;£ 



(l + £,7(^,^))^+^(0,0)^ 



p(0,0) = O 



G 



*(p(0,0),y(0,^),A,0;e) 



(l+£r7(0,</)))^+C.(0,<^)^ 



(5.35a) 



(5.35b) 



y(.9,i 
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by inserting rj{9, (p) from (|5.33p into and respectively. In [551 122] it is proved 
that the system of equations (|5.35p and (|5.28p has a locally unique solution (A,p, y), 
for \e\ sufficiently small, by considering the iteration 



2£2/3h(A*) - 



y 



(5.36a) 



(5.36b) 



;p('=+i)(0,0),l))=O, (5.36c) 

where 5X^''+^^ = A^'^+i) - A^^') and 

f(^-) {6, (p) ^ 2e2/3,(A*)p('=) {e, 4>) - a*{4>) ■ (p^^^^ {9, cb), y^"^ {9, 4>),\^''> , 
7^''^ {9, <t>) ^ {0, ct>) - G* f p(^) (0, 0), (0, 0), AC^-) , 9 



with starting values 

(Note that r/C") (6*, </>) and w('')(6', 0) are defined through and ((Oi)) respec- 

tively, using the values A^''^ p'^''\9,(j)) and y^^\9,(j)).) The key idea behind showing 
that these iterates remain bounded and then converge is to integrate (I5.36ap against 
(5/^ (^)^ after which the left-hand side becomes 

2e^(3^{X*) + £^r;(^-) (0, <^) + ^c^^) {9, 0)| p^'^+i) (0, 0), p^^^+i^ {9, cj,) 

Since (j5.32bp shows that the leading non-constant term in uj^^\9,(f)) is 0(£^), and 
(j5.32cp together with (|5.3ip shows that the leading non-constant term in ■q'^^^ (9, (j)) is 
O(e^), Assumption 15.51 ensures that the last expression is a definite quadratic term in 
(/)) for \e\ sufficiently small and this is sufficient for 28J to prove the following 

theorem. 

Theorem 5.4. Suppose F in p.ip has r > 5 continuous derivatives for (A, a;) 
in a neighbourhood of {X*,v*{9)). Then 3er > such that for \e\ < £r (|5.35p has a 
locally unique solution 

A*(£), p*(0,<^;£), y*{9,cj,-e) 

with p*(.,.;e), y*{.,.;£) having (r — 1) Lipchitz continuous derivatives. This means 
that both y and z, 'z ( through (|5.16p ) have this degree of smoothness, and so, through 
(j5.3p . do the invariant tori as manifolds. 

As in Theorem 14.41 in general ^ as r — * oo and we cannot expect analytic tori. 
In practice we seek an approximate solution of (j5.35p in the form 

p,,,,(0,0)=^p,,™e'(^''+"^) 

\£\ < L, \m\ < M 
y,,„(0,0)^5:y,,„e'(^«+™^) ' 
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with p_i _„i, y^i,^rn the conjugates of p^^m, Vi^rn- These functions must satisfy 

(5.37a) 



d_ 

d0 



Pl,m( 







Pofl = 0, 



(5.37b) 



(5.37c) 



for 1^1 < L, |m| < M and Va; e M"-^, with ryC^) (6*, 0) and w('=)(6',0) defined through 
(|5.33p and (|5.34p respectively, using the values A, Pl.m{0, (f>) and y^^ j^j{9, 0). As shown 
in 30J, the analogous iteration to (|5.36p also converges here and gives the following 
theorem. 

Theorem 5.5. Under the conditions of Theorem \5.4\ (|5.37p has a locally unique 
solution 



satisfying 



max ■ 



As in 



|Al,,(e)-A*(£)| 
llpf,M(->-;£) -Si„MP*(-,-;e)||L 

|yf,M 



(., .;£) - Si,My*(., •;e)||Li 
we comment on the implementation of this algorithm in 



\{\-S,.,,)p%,-e)\\H^ 
1(1 -Si.M)y*(.,.;£)|[Hi 



5.4. Numerical results. We consider a numerical example for which a group 
orbit structure leads to an interesting simplification of the general Neimark-S acker 
bifurcation equations: this is the Kuramoto-Sivashinsky equation in the form 



i<9 



^(x,t) + i-(.^(x,t)) 



(5.38) 



with u being both 27r-periodic and having zero mean in x [m 131] . We immediately 
obtain a finite-dimensional autonomous system by restricting to the Fourier approxi- 
mation 



l{x, i) « ^ uii{t)e 



ife uo{t) ^ 



and making use of the conjugacy condition leads to the complex system 
^ = -4DiM + A [DiM-iDLQ„(M)] mgC^, 



(5.39) 



(5.40) 



where is the Lx L diagonal matrix with entries 1, . . . , L and the quadratic function 
Q^:C^ ^ is defined by 



e-i 



(5.41) 
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(Thus (|5.4ip is our only discretisation error.) We describe below the sequence of com- 
putations which leads to Neimark-Sacker bifurcation for (|5.40p : this mimics some 
of the numerical results in [15', which should be referred to for further informa- 
tion. These computations exhibit our fundamental Crandall-Rabinowitz formulation 
in three different bifurcation situations. 

a) Bifurcation from the trivial solution. (|5.40p has the trivial stationary solution 
curve M = VA, which is stable for A < 4, and nontrivial stationary solutions bifurcate 
at 

Atb = 4£2 £=1,2,...,L. (5.42) 

These nontrivial stationary solutions are not isolated, since the autonomous nature 
of (|5.38p implies that if u{x,t) is a solution then so is u{x + a,t) Va G R. Conse- 
quently, in order to apply the Implicit Function Theorem and Newton's method, we 
must eliminate this multiplicity by either imposing a phase condition or a symmetry 
restriction. Since we are interested in the first bifurcation branch, i.e. X^y^ = 4 in 
(|5.42p . it is simplest to consider only stationary solutions of the form 

16 = is with s e M-^ 

for (|5.40p : this leads to the real system 

- 4Dis + A [Dls + DlQAs)] - 0, (5.43) 

where the quadratic function : is defined by 

For small we move onto the bifurcating curve of nontrivial stationary solutions 
by seeking solutions of (|5.43p in the Crandall-Rabinowitz formulation s = es, with 
amplitude condition si — 1. Hence, with starting values 

A(0)=4, s(°'=ei, 

the iteration in [6^ can be written 

A(fe+i) = 4 - er('^\ sf +i) = erf V(4£=^ - A(^'+')) e^2,...,L, 

where 

r« =AWg,(s«). 

b) Continuation of stationary solutions. Having moved away from the bifurca- 
tion point at Ajjj = 4, we can follow the branch of nontrivial stationary solutions 
by applying a standard continuation algorithm [T] to (|5.43D . This branch is always 
parametrisable by A, and so we can refer to solutions of (|5.43p by (A, s(A)) and the 
Jacobian matrix at solutions by 

Jf (A) = -4Di + A [Di + DiT.(s(A)) - DiH(s(A))] , (5.44) 

where in Matlab notation 



T,(s) = toeplitz([0 s{l:L~ 1)], [0 - s(l : L - 1)]) 
H(s) =hankel([s(2 : L) 0]). 
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L 


A* 

rw 


c* 






8 

16 


13.0038442196 
13.0038442196 


0.9990409957 
0.9990409957 


17.3973078781 
17.3973072209 


3.3475479311 
3.3475479124 



Fig. 5.1. Approximation of bifurcation points for the Kuramoto-Sivashinsky equation 



The eigenvalues of (|5.44p remain strictly in the left-half plane but this matrix, how- 
ever, only measures the effect of symmetric perturbations. To consider the effect of 
symmetry-breaking perturbations we must monitor the matrix 



J-(A) EE -4Di + A [Di + DiT.(s(A)) + DiH(s(A))] 



(5.45) 



which always has a null-vector Dis(A) because symmetry was imposed specifically 
to eliminate non-isolated stationary solutions. As A moves away from X*^ = 4, all 
the other L — 1 eigenvalues of (|5.45p remain at first strictly in the left-half plane 
but, as A approaches X*^ ~ 13, our zero eigenvalue becomes defective, with algebraic 
multiplicity two. At this value of A, we denote the null- vector of (|5.45p by e''^, with 
normalisation ||e™|| — 1, and the generalised eigenvector by cr™, with normalisation 
«^)2 + ||a™||2 = l, where 



c*.^Dls{X*.^ 



with 



.rw.T ^rw 



0. 



As part of our continuation algorithm, we can monitor the real part of the eigenvalues 
of (|5.45p and detect a crossing of the imaginary axis: a simple secant iteration then 
accurately determines the value of A at which bifurcation occurs and this is displayed 
in Figure lSTTI We can also check that the crossing is transversal, by using a simple 2nd- 
order centered finite difference (with step h) to obtain the following approximations 
to the critical eigenvalue derivative. 



h 


0.1 


0.01 


0.001 


Eigenvalue speed 


6.127497 


6.127414 


6.127418 



c) Bifurcation to rotating waves. This loss of stability is associated with the 
creation of a special type of periodic orbit called a rotating wave. It is a solution 
of with ([Og]) having the form 



where the unknown wave-speed c G M plays the role of "frequency" . The important 
practical point is that these rotating waves are as easy to compute as stationary 
solutions, since under the moving frame 



they satisfy 



ct 



d 



0, 



(5.46) 



where now v is 27r-periodic and has zero mean in ^. Hence, instead of (|5.39p . we use 



(5.47) 
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and arrive at the complex system 

- 4D> + A [Dlv - iDLQciv)] - icDlv = v eC^, 



(5.48) 



which is the analogue of (|5.40p . We can then move onto the curve of rotating waves 
by seeking a solution of (|5.48p in the Crandall-Rabinowitz formulation 

c = ec and v = is(A) + e [v^ + ivi] Vj^ , e 

for small |e|. Just as for ordinary Hopf bifurcation, we must complement (|5.48p with 
amplitude and phase conditions, and these are 

(<T™f i;, + <^c=l and (e™)^7;, = 0. 

Thus, splitting (j5.48p into real and imaginary parts, our analogue of the Hopf bifur- 
cation iteration in section Sj^l is 



' R 



1 



and Jf(AW)«^'+'^ =£r 



where 



) = Dis(A(^)) and = [r^)' [X^Jv^^^ + c«Dis'(A:^) 

r« ^ c«D,^('='+AWD^Re{Q,(^W+it;('^))}, 
with starting values 



A(°) = a; 



s(o) 



,(0) 



d) Continuation of rotating waves. Having moved away from this pscudo-Hopf 
bifurcation point A*^, we can follow the branch of rotating waves by applying a 
standard continuation algorithm [1] to (|5.48p . This branch is parametrisable by A and 
so, splitting V into real and imaginary parts, we can refer to the solutions of (|5.48p by 
(A, c(A), t;n(A) + it>i(A)). The analogue of Floquet exponents for the rotating waves 
are the eigenvalues of 



J™ (A) 



J™ (A) J™(A)^ 



i,2Lx2L 



(5.49) 



where 



JrT(A) 

J™ (A) 
J™ (A) 



1 



c(A) 
A 

~x 
1 



{-4Di + ADi + XDl [TMX]) + H{v,iX))]} 
Dl [T.(-Uh(A))-H(-(;,(A))] + Dl 

{-4Di + ADi + XDl [T.(^:(A)) - H(t;,(A))]} 
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and in Matlab notation 

T,(vr) = toeplitz([0 v4l : L - 1)]). 
As with periodic orbits, one of these is always zero since 

e™(A) 









e7(A) 







is a null- vector because of the autonomous nature of (|5.46p . Thus (|5.48[) must be 
complemented by a phase condition 

(e™(Aprcv))^t^K + (er(Aprcv))^^, - 0, 

where e™(Aprcv) is obtained from the solution at the previous value of A. Apart from 
this, all the other 2i — 1 eigenvalues of (|5.49[) lie strictly in the left-half plane until 
A approaches A*^ ~ 17.4, when a complex-conjugate pair ±ia;*g cross the imaginary 
axis with the complex eigenvector satisfying 



for 



As part of our continuation algorithm, we can monitor the real part of the eigenvalues 
of (|5.49p and detect a crossing of the imaginary axis: a simple secant iteration then 
accurately determines the value of A at which bifurcation occurs, with 



" ns' 






. I 







denoting the null-vector there. The variation with A (in the upper-half of the complex 
plane) of this critical complex-conjugate eigenvalue is shown in Figure 15. 2|, which 
may be compared with Figure 4.2 in |15j . while numerical values for Neimark-Sacker 
bifurcation are displayed in Figure [STTl (For L — 32, all results agreed to 10 decimal 
places.) As with bifurcation to rotating waves, we can also check that the crossing is 
transversal by calculating the following approximations to the real part of the critical 
eigenvalue derivative. 



h 


0.1 


0.01 


0.001 


Eigenvalue speed 


0.0966405 


0.0966218 


0.0966216 



e) Bifurcation to invariant tori. We seek invariant tori in the Crandall-Rabinowitz 
formulation 



M 



z-e, 



Vm 

Zi^m V^, TO 



for small so that T, together with tj, 77 G R, solves the finite-dimensional restriction 
of 



1 r d^T ^ 
-4-— r - A 



c(A) 



\- - — T ) 

Q^2 ^ 2Q^\ ) 



,dT dT ^ 



Here (A, c(A), Vh(A) -I- ii;i(A)) satisfy (|5.48p and, if we denote the (/)-modes of z by 
2; (to) G C^, then by conjugacy we need only solve for z{Q) = Zb{Q) + i2;i(0) with 
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3.75 r 



3.7 ■ 



3.65 



3.6 ■ 



3.55 



3.5 ■ 



3.45 



3.4 ■ 



3.35 



3.3 ■ 



-A,= 14 



X=18 



3.251- 
-O.f 



-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 

Fig. 5.2. Movement in C for critical eigenvalue of rotating waves 



ZaiO), Zi{0) G and z+{m), z_{m) e for m = 1, . . . , M. (Note that z+ contains 
the ^- modes for £ > 1 and 2;_ for ^ < — 1.) Hence, dividmg through by e, we may 
write the finite-dimensional restriction as 



J™(A) J;7(A) -DLV^iX) 



ZniO) 
V 



MO) 

r,(0) 



for m = and 



for m > 1, where 



[J™ (A) - iujml] 



jr(A) 



z+(m)' 
z_ (m) 



r_|_(m) 
r (m) 



J™ (A) J™(A)- 
J™ (A) J™ (A) 



is the complex version of (I5.49P defined by 
1 



J™ (A) 
J™ (A) 



c(A) 



{-4Di + ADi - iADiT.(t;,(A) - ivAX))} - iDl 



-i-^DiHK(A)+i-u,(A)) 
c(A) 

J™ (A) = i^DiHK(A) - iv,{X)) 



J™ (A) 



c(A) 
1 



{-4Di ADi -t- iADiT,K(A) + iv,{X))} + iD^. 
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The right-hand sides are defined by 



and 



r-a(O) 
r,(0) 



r+(m) 
r_(m) 



A 



c(A) 



-DLq_{m) 



+ V 



+ it] 



DiZR(O) 



-DiZ_(m) 



with qj^(0),qj(0) e and q^{m),q_{m) G being derived from the quadratic 
term 



is 



m=-M 



qo,m = Vm 



in the same way as their analogues for z. We also have the phase condition 

(e-f ZR(0) + (errz:(0)=0 

and, since 

we can define our amplitude and phase conditions using 





1 ,ns" 




0" 




: _ns 
-ICTj 








where a € C is chosen so that 



icr;- 
icr" 



iKir +iicT^ii' = i. 



Thus our Neimark-Sacker bifurcation iteration is 



■Ja^CA^'^)) JaTCA^^)) 

J™(AW) J™(AW) 



DLt;:(AW) 
-Dii;R(AW) 
















= = e 














for m = 0; 



irw 



(A«) 
J™(AW) 



for m = 1; and 



(fc)l 



irw 



(AW) 
(AW)-ia;W| 







J™(A('=))-imw 



(fe)| 



'z^^+'\m) 




r''^^ (to) 


z^}+'\m)_ 


= e 


r_ (to) 
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for m > 2. Here 











~ c(AW) 





r\ (m) 
r_ (m) 



DLqi'\m) 
-DLqi^\m) 



and we note that, for m = 1, our two extra real unknowns uj, A are compensated by 
one extra complex condition. Our starting values are 

A(°) = a:,, ^(°) = ^(«) = 0, - - -"^ 

with all other components of z^"-* zero. In Figure [5731 we display the numerical results 



V 




A 


c(A) 


-0.1652838896 


3.1325935189 


17.3749646320 


16.7931284737 



11^(0)1! 


0.25 


11^(1)11 


1.00 


11^(2)11 


0.10 



IN(3)|| 


0.11 X 10"^ 


Mm 


0.13 X 10-^ 


M5)\\ 


0.15 X 10-^ 



11^(6)11 


0.18 X 10"* 


11^(7)11 


0.21 X 10"^ 


11^(8)11 


0.25 X 10"^ 



Fig. 5.3. Decay of Fourier 4>-modes for invariant torus at e = 0.5 



for L = 16 and M = 8, in particular verifying the decay of the size of the (/)-modes 
for = [z^,z'^) . 

f) Concluding remarks. Finally, we emphasise the simplification in the Neimark- 
Sacker algorithm that the group orbit structure of the Kuramoto-Sivashinsky equation 
allows. Just as the rotating waves are really periodic orbits that can be calculated 
as simply as stationary solutions, so the invariant tori can be calculated as simply as 
periodic orbits: i.e. there is only one explicit independent periodic variable and thus 
no resonance can occur. This means that we can utilise the simple parametrisation for 
the tori in H5.ll (as above with constant t] and ui) for an arbitrary number of i^-modes, 
rather than being limited to M = 3 by Assumption 15.41 

6. Conclusion. In [JT] we stated that the fundamental idea behind the present 
paper is to use the approach in \28\j . . .to develop a practical computational algorithm 
for Neimark- Backer bifurcation. We claim to have achieved this goal, but the final 
implementation of the algorithms in i J4.4l and ij5.3l will be explored elsewhere. The 
two main reasons for this are the length of the present paper and the belief that these 
practical questions are best-suited to a separate paper. We emphasise, however, the 
two key points that an efficient algorithm must address. 

a) The extraction of normal form information from the simple low-order Fourier 
approximations in ij4.2l and i j5.2l which then allows us to introduce the es- 
sential, but more complex, parametrisations in H4.4I and >^5.3I 
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b) Our final iterations in (|4.35p , (|4.37p , (|5.36p and (|5.37p necessarily rely on the 
solution of linear variable- coefficient differential equations. This raises the 
question of computational efficiency since, throughout this paper, we have 
utilised the mode-decoupling property for Fourier approximations of constant- 
coefficient systems. Our solution to this problem is to make use of the precise 
structure of the variable-coefficient equations in order to pre-condition them 
by suitable constant-coefficient operators [H HJ [32] . 

Finally, we remark on several other points which, for the sake of simplicity, were 

omitted earlier. 

• In [2] and !j5] we assumed that the basic periodic orbit v*{9) was known 
exactly. In practice, of course, we would have a sufficiently accurate Fourier 
approximation, as in iJ4.3l and i J5.4l 

• In [21 and fj5]we assumed that — for the Floquet theory described in ^31 
The case n_ > introduces no practical difficulties, whether these eigenvalues 
occur in E* or E*. In both cases, the strategy in [22] can be followed. 

• We have avoided any discussion of aliasing, numerical quadrature and the 
FFT for our Fourier spectral methods [21 [4] , by implicitly assuming that all 
integration was performed exactly. The only practical difference is that some 
of our errors in Theorems 12. 4[ 14.31 and 15.31 may be 0{e) rather than 0{£^). 
This is still sufficient for our purposes, but may be avoided if desired: such 
questions will be addressed in the future paper mentioned above. 

• We have merely stated the smooth invariant subspace decompositions re- 
quired in (|2.2p , (|4.3p and (|5.2p . Further information may be found in [7] . 

REFERENCES 

[1] E.L. Allgower and K. Georg, Numerical Continuation Methods: An Introduction, Springer, 
1990. 

[2] D.J. Allwright, Harmonic balance and the Hopf bifurcation, Math. Proc. Camb. Phil. Soc, 

82 (1977), pp. 453-467. 
[3] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Dover, second ed., 2001. 
[4] C. Canuto, M.Y. Hussaini, Quarteroni a., and Zang T.A., Spectral Methods in Fluid 

Dynamics, Springer, 2003. 
[5] P. Coullet and E. Spiegel, Amplitude equations for systems with competing instabilities, 

SIAM J. AppL Math., 43 (1983), pp. 776-821. 
[6] M.G. Crandall and P.H. Rabinowitz, Bifurcation from simple eigenvalues, J. Funct. Anal., 

8 (1971), pp. 321-340. 

[7] J.W. Demmel, L. Dieci, and M.J. Friedman, Computing connecting orbits via an improved 

algorithm for continuing invariant subspaces, SIAM J. Sci. Comp., 22 (2000), pp. 81—94. 
[8] E. Doedel and L.S. TuCKERMAN, eds.. Numerical Methods for Bifurcation Problems and 
Large-Scale Dynamical Systems, vol. 119 of The IMA volumes in mathematics and its 
applications. Springer, 2000. 
[9] J. GuCKENHEIMER and p. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurca- 
tion of Vector Fields, Springer, 1983. 

[10] B. Hassard, N. Kazarinoff, and Y.-H. Wan, Theory and Applications of Hopf Bifurcation, 
Cambridge University Press, 1981. 

[11] B. Hassard and Y.-H. Wan, Bifurcation formulae derived from center manifold theory, J. 
Math. Anal. AppL, 63 (1978), pp. 297-312. 

[12] G. lOOSS and M. Adelmeyer, Topics in Bifurcation Theory and Applications, no. 3 in Ad- 
vanced Series in Nonlinear Dynamics, World Scientific, 1992. 

[13] G. lOOSS, A. Arneodo, P. Coullet, and C. Tresser, Simple computation of bifurcating 
invariant circles for mappings, in Dynamical Systems and Turbulence, D. Rand and L.-S. 
Young, eds., no. 898 in Lecture Notes in Mathematics, Springer, 1981, pp. 192-211. 

[14] G. lOOSS AND D.D. Joseph, Elementary Stability and Bifurcation Theory, Springer, second ed., 
1990. 



NEIMARK-SACKER BIFURCATION 



43 



[15] I.G. Kevrekidis, B. Nicolaenko, and J.C. Scovel, Back in the saddle again: a computer 
assisted study of the KuramotoSivashinsky equation, SIAM J. Appl. Math., 50 (1990), 
pp. 760-790. 

[16] B. KrauSKOPF and H.M. Osinga, Investigating torus bifurcations in the forced Van der Pol 
oscillator, in Doedel and Tuckerman [8], pp. 199-208. 

[17] Y. A. KuzNETSOV, Elements of Applied Bifurcation Theory, Springer, third ed., 2004. 

[18] O.E. Lanford, Bifurcation of periodic solutions into invariant tori: the work of Ruelle and 
Takens, in Nonhnear Problems in the Physical Sciences and Biology, I. Stakgold, D.D. 
Joseph, and D.H. Sattinger, eds., no. 322 in Lecture Notes in Mathematics, Battelle Sum- 
mer Institute, Seattle, Springer, 1973, pp. 159-192. 

[19] J. Marsden and M. McCracken, Hopf Bifurcation and its Applications, Springer, 1976. 

[20] G. Moore, Geometric methods for computing invariant manifolds. Applied Numerical Math- 
ematics, 17 (1995), pp. 319-331. 



[21] , Computation and parametrisation of invariant curves and tori, SIAM J. Num. Anal., 

33 (1996), pp. 2333-2359. 

[22] , Floquet theory as a computational tool, SIAM J. Num. Anal., 42 (2005), pp. 2522-2568. 

[23] J. MOSER, On the theory of quasiperiodic motions, SIAM Review, 8 (1966), pp. 145-172. 

[24] , Stable and Random Motions in Dynamical Systems, no. 1 in Princeton Landmarks in 



Mathematics, Princeton University Press, 2001. 

[25] J. MuRDOCK, Normal Forms and Unfoldings for Local Dynamical Systems, Springer, 2003. 

[26] V. Reichelt, Computing invariant tori and circles in dynamical systems, in Doedel and Tuck- 
erman [8], pp. 407-437. 

[27] D. Ruelle and F. Takens, On the nature of turbulence. Comm. Math. Phys., 20 (1971), 
pp. 167-192. 

[28] R.J. Sacker, On invariant surfaces and bifurcation of periodic solutions of ordinary differential 

equations, IMM-NYU 333, Courant Inst. Math. Sci., New York Univ., 1964. Available at 

http : //www-rcf .use . edu/'rsacker/pubs/IMM-NYU.pdf . 
[29] , A new approach to the perturbation theory of invariant surfaces. Comm. Pure and 

Appl. Math., 18 (1965), pp. 717-732. 
[30] A.M. Samoilenko, Elements of the Theory of Multi- Frequency Oscillations, Isdatelstvo Nauka. 

Moskva, 1987. English translation: Kluwer Academic Publishers, Dordrecht et al., 1991. 
[31] R. Temam, Infinite- Dimensional Dynamical Systems in Mechanics and Physics, no. 68 in 

Applied Mathematical Sciences, Springer, 1988. 
[32] H.A. VAN der VorST, Iterative Krylov Methods for Large Linear Systems, no. 13 in Cambridge 

Monographs on Applied and Computational Mathematics, Cambridge University Press, 

2003. 

[33] Y.-H. Wan, Bifurcation into invariant tori at points of resonance. Arch. Rat. Mech. Anal., 68 
(1978), pp. 343-357. 

[34] , Computations of the stability condition for the Hopf bifurcation of diffeomorphisms on 

R2, SIAM J. Appl. Math., 34 (1978), pp. 167-175. 



